program poisson

    use rhofunc, only: rho
    use command_line, only: get_m_value_from_CLI
    implicit none

    integer, parameter :: dp = kind(0.d0)
    integer :: i, j, iter, M
    real(dp), parameter :: epsilon0 = 8.85E-12_dp, target = 1E-6_dp, a = 0.01_dp, a2 = a*a, &
                           const = a2/4.0_dp/epsilon0
    real(dp) :: delta, b, e
    real(dp), allocatable :: phiprime(:, :), phi(:, :), rhoarr(:, :)

    M = get_m_value_from_CLI()
    allocate (phiprime(M, M), phi(M, M), rhoarr(M, M))
    phiprime = 0.0_dp
    phi = 0.0_dp

    delta = 1.0_dp
    iter = 0
    call cpu_time(b)

    do i = 1, M
        do j = 1, M
            rhoarr(i, j) = rho(i*a, j*a)
        end do
    end do

    do while (delta > target)
        iter = iter + 1

        do i = 2, M - 1
            do j = 2, M - 1
                phiprime(i, j) = (phi(i + 1, j) + phi(i - 1, j) + phi(i, j + 1) + phi(i, j - 1))/4.0_dp &
                                 + const*rhoarr(i, j)
            end do
        end do

        delta = maxval(abs(phiprime - phi))
        phi = phiprime

    end do

    call cpu_time(e)
    print '(a,g0.4,a,g0.4)', "time elapsed: ", e - b, " s, iteration: ", iter

end program
